Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 28 April 2010 (MN JAT^ style file v2.2) 



The dynamics of pulsar glitches: Contrasting 
phenomenology with numerical evolutions 

o "■ 

^ T. Sidery^, A. Passamonti^ and N. Andersson^ 

FENS, Sabanci University, Orhanli, 34956 Istanbul, Turkey 
5-H ^ School of Mathematics, University of Southampton, Southampton SOU IBJ, United Kingdom 



o 



m 
d 



o 



28 April 2010 



(N 

W ■ ABSTRACT 

I I I , In this paper we consider a simple two-fluid model for pulsar glitches. We derive the 

basic equations that govern the spin evolution of the system from two-fluid hydrody- 
^ ^, namics, accounting for the vortex mediated mutual friction force that determines the 

glitch rise. This leads to a simple "bulk" model that can be used to describe the main 
properties of a glitch event resulting from vortex unpinning. In order to model the 
^ ' long term relaxation following the glitch our model would require additional assump- 

, tions regarding the repinning of vortices, an issue that we only touch upon briefly. 

Instead, we focus on comparing the phenomenological model to results obtained from 
, time-evolutions of the linearised two-fluid equations, i.e. a "hydrodynamic" model for 

^ ' glitches. This allows us to study, for the first time, dynamics that was "averaged" in 

OO , the bulk model, i.e. consider the various neutron star oscillation modes that are excited 

during a glitch. The hydro-results are of some relevance for efforts to detect gravita- 
tional waves from glitching pulsars, although the conclusions drawn from our rather 
simple model are pessimistic as far as the detectability of these events is concerned. 
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1 INTRODUCTION 

Neutron stars provide excellent testbeds for extreme physics theory. Since their core density reaches beyond anything we can 
produce in the laboratory, observat ions of such compact remnants may provide unique constraints on models for supranuclear 



physics (jLattimer fc Prakashll2004l ). In order to infer useful information from astrophysical observations we need to construct 
realistic models with the power to predict the evolution of individual systems (at least to some extent). This is a serious 
challenge for the modelling community. Even a moderately reasonable neutron star model should account for the presence of 
different exotic states of matter. From nuclear physics and Bardeen-Cooper-Schrieffer (BCS) theory we expect the outer neu- 
tron star core to consist of superfluid neutrons, superconducting protons, free electrons and muons. Deeper into the core more 
exotic phases of matter, like superfluid hyperons and/or deconfined quarks exhibiting colour-flavour-locked superconductivity, 
are likely to be present. Meanwhile, a relatively thin (1 km or so) crust surrounds the fluid core. In the crust, matter changes 
from a soup of nucleons in the interior to an elastic nuclear lattice of heavy iron near the surface. In the inner part of the 
crust (beyond neutron drip) free neutrons are expected to be superfluid. 

An important question concerns whether differences in internal structure can be deduced from observations. This is 
obviously problematic since most of the data collected for these stars originate from the star's surface, atmosphere or mag- 
netosphere. Having said that, there are phenomena that involve bulk dy namics and which should depend on the internal 
composition. Examples of such phenomena are the radio pulsar gli tches ( Lyne et al.| 200C ) and the quasi-periodic oscilla- 



tions observed in the X-rays from magnetar flares, see for example Watts fc Strohmaveil ( 2007 ). The long relaxation time 
associated with g litches is seen as indirect evidence for neutron star superfluidity ijAnderson fc Itoh 19751 : Ruderman 19761 : 



Alpar et al.l Il98ll 'l. and recent considerations of the crust oscillations which are thought to be the origin of the observe d 



magnetar oscillations suggest that these may be affected by the presence of a crust superfluid as well ( Andersson et al. 20091 ) 
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In this paper we discuss basic models for pulsar glitches. We focus on the glitch event itself rather than the subsequent 
long term relaxation. We consider the implications of the standard two-fluid model from two different points of view. First 
of all, we derive the simple equations that govern the bulk evolution of the system from two-fluid hydrodynamics, account- 
ing for the vortex mediated mutual friction force. This l eads to a basic mode l that can be used to describe a glitch rise 
resul t ing from global vortex unpin ning, be it in th e crus t (Li nk fc Epstein 19961 : Melatos et al. 20081 : Warszawski fc MelatosI 
20081 : iMelatos fc Warszawskilboogl l or in the core |Linkll2003l lT The results demonstrate that the model requir es additional 



assum ptions regarding the repinning of vortices in order to model the long-term evolution. This is as expected (|Alpar et al 



1984 ). Secondly, we use time-evolution of the linearised two-fluid equations as a "hydrodynamic" model for the glitch event. 



This allows us to consider dynamics that was "averaged" in the bulk model. In particular, we consider the various neutron 
star oscillation modes that are excited during the glitch. In principle, the obtained results could be of relevance for efforts to 
detect gravitational waves from glitching pulsars. Having said that, our estimates suggest that the gravitational-wave signals 
associated with the impulsive events that we consider here are very weak. 



2 BULK PROPERTIES 



We want to consider the simplest viable model for large pulsar glitches. The angular momentum of any superfluid component 
is determined by the density and configuration of vortices threading the fiuid. If the vortices are fixed (pinned), there can be no 
angular momentum exchange between the superfiuid and other components in the star. Assuming a scenario of catastrophic 
unpinning it is straightforward to formulate a simple glitch model. One can simply assume that the system has two components, 
with different moments of inertia and spin-rates. Adding the assumption that vortex pinning allows a lag in the rotation 
between the observed component (e.g. the crust) and the other component (the interior superfluid) and that the pinning 
breaks once a critical lag is reached, one arrives at a phenomenological glitch model. In order to effect the rapid transfer 
of angular momentum that leads to the observed spin-change in the crust, one would typically assume that the rotation 
lag relaxes on some timescale. The s t andard model assumes that this rela xation is due to the mutual friction acting on the 
superfluid vortices ( Alpar et al. 1984 : Mendel] 1991 : Andersson et al. 20061 '). 

While such a phenomenological model is useful, its scope is limited. Ultimately, a detailed description will require a 
hydrodynamics anal ysis. In particular , if w e want to understand the reason f or the sudden vortex unpinning (likely due 



to an instability, see Andersson et al. ( 2003t ): Glampedakis fc Andersson (200^) for recent ideas) and possible neutron star 



oscillations excited by the event. Developing a detailed hydrodynamics model is a severe challenge given the many uncertainties 
in the relevant physics, but we can make some progress by making contact between the general two-fiuid hydrodynamics 
framework and the phenomenological bulk dynamics description. As a suitable starting point, we will show how the bulk 
model can be obtained from hydrodynamics. This is instructive since it provides insight into the validity of the model, and 
also gives us a better idea of the origin of the different parameters (like the global spin- up time). Moreover, we will be able 
to make direct comparisons with hydrodynamics results. 



2.1 Single fluid 

It is useful to begin by outlining the analysis of a single fluid body. In that case, we have a velocity field Vi which evolves 
according to the Euler equations; 

(^^ + «'v,J«. + v.(A + 0) = o, (1) 

where 4> is the gravitational potential and jl = fj,/m is the chemical potential divided by the particle mass. In addition, we 
have the continuity equation 

|+V.(p.^)=0, (2) 

and the Poisson equation 

V^<;/> = 47rGp. (3) 

Note that we use a coordinate basis to represent vector quantities throughout this paper. In other words, we distinguish 
between co- and contravariant quantities, using the flat space metric gij to relate them. That is, we have Vi = QijV-' . We also 



^ In reality, our model is somewliat unrealistic for both crust and core superfluids. In the former case we have not accounted for the 
crust elasticity, while in the latter case we are ignoring the expected interaction between neutron vortices and proton fiuxtubes. These 
effects may have a decisive impact on glitch dynamics. Nevertheless, the present work is state-of-the-art in this area, and we expect to 
add the relevant features to the hydrodynamical model in future work. 
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make use of the Einstein summation convention for repeated indices. Finally, we have assumed that the internal energy E 
depends only on the number density n, i.e. that the fluid is a barotrope. Then 

and the pressure P of the fluid is defined (in the usual way) by 

dP = ndfi . (5) 

Viscosity terms have been omitted in anticipation of applying the solid body approximation, which we do now. Assuming 
uniform rotation the fluid velocity can be written 

Vi = tijk^^ = Vlzuipi , (6) 

where we take f2 = vj is the cylindrical distance from the rotation (z) axis and (fii is a unit vector in the direction of the 

flow. By assuming axi-symmetry it follows that 

v'VjV.^O, (7) 

vi\/,ifi + (k)^0. (8) 

We will now derive the conservation of angular momentum and energy from the above equations. Contracting ((T} with pvi 
and integrating over a fixed volume V gives 

This shows that the kinetic energy is conserved. Adding the assumption that the fluid is in solid body rotation we see that 

v"" = QjQ'' (^Slx^ -x'xk^ , (10) 
where x^ = QijX^x-' = x-'xj. From this we define the moment of inertia as 

= j pi^Slx" -x'xi)dV . (11) 
We can use (|10p and (lllfl to rewrite the change in kinetic energy equation ^ as 

dt 2dt\^ ) ^ ' 

Let us now consider the z component of the angular momentum. Assuming that the chemical and gravitational potential 
only depend on the (spherical) radial position, i.e. assuming slow rotation, then 

eijkX^V'' {jl + (j)) — Q for i — z. (13) 
Contracting £i with peijkX-' and integrating gives 

dJ, 
dt 



^ = / petjkX^S'^dV 



d_ 

dt 



I pQj (^Six'^ ~ Xix'^ dV (I'^j) =0 for i = z. (14) 



We see that, for cylindrical polar coordinates with Q,^ — (0, 0, Q) and — I we get the standard results 
Both these quantities are conserved. 



E=^m'^ and J" = IQ . (15) 



2.2 Two-fluid model 



We now consider a two constituent stellar model. Our particular interest concerns two effects, the entrainment and the mutual 
friction between the components. In order to simplify the initial analysis, we first ignore the mu tual friction . 

Our formulation for multi-fluid hydrodynamics in Newtonian gravity derives from the work by Prixl ( 20o3 l . see also Andersson fc Comer 



(2001 



The analysis is based on a variational principle, where the action is minimised by varying the fluid flow lines. Key in 



this analysis is the allowance that the internal energy may depend on the velocity difl^erence between the two constituents, 
wf^ = vf — vf, such that 



dE 



dE 



dUn + 



dE 



drip + 



dE 
dwip 



d™np = fJ-ndrin + pipdup + adw^p . 



(16) 



Here, the constituent indices n and p denote the neutron and "proton" (incorporating also electrons and muons in the usual 
way) components, respectively. An important consequence is that the conjugate momentum density for each constituent is 
modified so that 



px + Ex («f - -uf)] , 



(17) 
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where x and y is either n or p, with x 7^ y. In this relation the entrainment is represented by the parameter Ex- This is 
a non-dissipative effect, in the neutron star case due to the strong interaction between neutrons and pr otons, which lead s 
to the momentum n o longer being align e d wit h the individual compone nts velocity. We refer the read to iPrix et alj (|2004l ): 
Carter et all l|2005h : lGusakov fc Haense] |2005l ): lchamel fc Carter] |2006l ) for discussions of the role of entrainment in neutron 



star dynamics. Note that, 

ExPx = £yPy = 2a . 

We now have a set of Euler equations for each constituent. These equations can be written 

d 



dt 



0. 



(18) 



(19) 



In addition, we have one continuity equation for each component and the Poisson equation for is now sourced by the total 
mass density p = pn + Pp. Following the analysis in the previous section, we want to derive the equations that represent the 
global conservation of energy and angular momentum. In doing this we will, for simplicity, assume that the velocity fields of 
the constituents are those of rotating solid bodies. We also assume that the two components rotate around the same axis, i.e. 
we ignore any precessional motion. This means that we can write 



tlx ~ Q,^zoLp^ and Wyy^ = — vl^ = (S7y — 57x) zotp^ 



(20) 



As we are assuming solid body rotation. Six is not a function of position. Moreover, the axial symmetry of the system implies 
that 



wiV,(<A + /ix) = 0. 



(21) 

,,vv^ + /ix) = 0. (22) 

We continue to follow the method used in the single fluid case, contract £f with pxt'f and integrate over a fixed volume to 
find the global change in energy of each constituent. This leads to 

We obtain the final result for the total change in energy by adding the expressions for the two components, 



aSx 
dt 



p^vlEldV = i 



dV = 0. 



(23) 



dE 
'dt 



dE^ . dE^ 



dt 



+ 



dt 



ld_ 

2 dt 



4awi^\ dV} = 0. 



(24) 



This result shows how the entrainment affects the conserved energy, c f. Carter &: Chamell (1200J) for a detailed discussion, and 
agrees perfectly with the modified "kinetic energy" used bv lMendelll (|l99lf ). 

In the particular case of (aligned) solid-body rotation, with flf aligned with the z-axis, we have 

dE d 



dt dt 

Here we have defined the constituent moment of inertia as 



i/nfin [fin + En (flp - On)] + ^/p^p [fip + Ep (fin 



flp 



0. 



dV , 

and Jx = /xz. Similarly, we can calculate the total change in angular momentum. To do this we note that 



tijkX'' v!^\/''vf = eijkX'' x'^D,'^ — eijfe3;"'Slx3^;f2x = for 



(25) 
(26) 
(27) 



as flf is parallel to the z-axis. Contracting £f with px^ijkX-' and integrating over the volume V we arrive at 

dJl 
dt 



= / p^eijkX^ S'^dV 



-ex. 



dvt 



dt 



dV = for 



This can be rewritten as 



^ = ^{^xH«.^ + e.(f^I-«n]}=0 for 



(28) 



(29) 



from which it follows that the total change in angular momentum is given by 

^ = ^(/nfln+Jpflp) =0. 

Hence, the total angular momentum is conserved. This is obviously not surprising. The non-trivial result concerns how the 
entrainment affects the evolution of the individual components. This will be important later. 



(30) 
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3 A SIMPLE SPIN-DOWN MODEL 

We are not yet in a position where we can model glitches. To do this, even in the most basic fashion, we need to account for 
the coupling due to mutual friction. However, before we discuss that problem it is interesting to consider how the conservation 
equations that we have obtained can be used to model the rotational evolution of the system. The main purpose of doing this 
is to understand how the entrainment enters the problem. As we will see, the result can be quite surprising. 

Equation (|29|) shows that the angular momentum of each constituent is conserved. In a neutron star we would expect the 
protons to be locked to the magnetic field and the crust. Hence, they should be spun down due to a magnetic torque J^^™. We 
can include this torque in the equations by breaking the conservation of the proton angular momentum, so that 

jf = - ir , (31) 

(from now on we will often represent time derivatives by dots). To be specific, we assume that the magnetic torque is related 
to the angular velocity of the protons by 

tTcni = vA/pilp . (32) 

This is in accord with the standard magnetic dipole model. We will also assume that the constituent moments of inertia are 
constant. Noting that /nEn = Ip£p and defining e = Ep equation (|29|) leads to 



Defining / = Ip/In equation (|33p gives 
Substituting this into equation (|34p we get 



1-eJ 

This is a separable equation so we can integrate to get 



IMn + IpE [Qp - fin) = , (33) 

Ipilp + IpE {tin - tip) = -AIp^l . (34) 



tin = - ^ '^^ - tip . (35) 



-Anl . (36) 



/■"p d£ip^_A r*^^ 

J do ^ Jto ' 



where Slo is r^p at time to. Setting to ~ we find the solution 



(37) 



!^p = fio (^1 + ^^-p-J • (38) 
As we expect the evolution to be slow, the second term in the bracket is small and we can expand to get 

Hp « Ho f 1 - ^) . (39) 



I 

From this we can find the characteristic evolution timescale r for the crust (the protons). Ignoring entrainment we have 

^ = :4pr 

A can be calculated from the standard magnetic dipole model ( Shapiro fc Teukolskv 19831 ). Modifying the result so that the 
torque acts on the proton fiuid rather than the whole star we find 

where Bp is the strength of the magnetic dipole with axis at an angle 6 to the rotation axis, R is the radius of the star and c 
is the speed of light. As 7p ^ In « / (typically), we can use 

/p.^/.^^, (42) 

In In 5 

where M is the mass of the star. Using typical parameters B = lO^^G, R — 10® cm, Ip/In ~ 0.05, Q.o = 27r/0.1 s^^ and 
M = 1.4M0 and setting = 7r/2 we find r « 7 x lO" years. 

Let us now consider the evolution of the, unseen, neutron component. From equation (|35l) we have 



fin - nl = ^ (fip - fio) , (43) 

I — el 



where Q.^ — Q.n{t — 0). Substituting equation (|39p into this we get 



nn~nl + - (44) 

1 - eJ T 
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This relation allows us to estimate how long it takes for a rotational lag to develop between the two constituents. Assuming 
that the two components rotate together at time t = then $7^ = fio. The rotational lag then evolves according to 



An = fin 



el 



= 1 + 



eI r 
ei 



T 



I -el 



T 



el 



(45) 



AD. 



t 



l-eir 

It has been arg ued (|Lvne et all boool ) that a rotational lag of AD/Qp 
glitches. From 1)46^ we see 



i « 10- 

r 



f ~ 7 years . 



(46) 

10~* is needed in order to "explain" Vela sized 



(47) 



Hence, this simple model is consistent with large glitches occuring once in a few years in a typical young pulsar. 

Finally, let us consider the entrainment coupling in more detail. In general, the evolution of the rotation of the proton 
and neutron fluids is given by p9p and (|44|l . respectively. If we focus on a sufficiently short evolutionary timescale, then we 
can assume that Jem is approximately constant. From equations p3|l - (|35|l we find 



n„ = -— 1- 



eln 



In - eh 



J CI 



(48) 



This is an interesting, and perhaps surprising result. It appears that, even though a spin down torque acts on the protons 
their rotational velocity may increase. This happens when 



In ^ ^ In 
In + Ip Ip 



(49) 



Woul d this happen for realistic parameter values? Setting e — 0.05 and Ip/In ~ 0.1 , typical values in the outer neutron star 
core (jChamellbooi ), we find that the crust should spin down. However, if we consider the neutron fluid we find the condition 
for spin up is 

< e < • (50) 

For the expected values of e and Ip/In we see that the neutron fiuid spins up. This is somewhat counter-intuitive, but does not 
violate any fundamental principles. In fact, it is easy to see how this effect can result from an exchange of angular momentum 
in a coupled system. Unfortunately, since we do not observe the neutron component directly, the result does not help us 
constrain the entrainment parameter. It is just an example of the drastic effect that entrainment can have on the evolution 
of a system. 



4 MODELLING A GLITCH 



The two-component model that we have described cannot be used to model glitch events unless we add more physics to it. 
Key to the problem is the motion of the superfluid vortices. There are two aspects to this problem; We need to account 
for the frictio n that arises due to the presence of the vortices and the associated dissipative coupling between the two fluids 
in the model ( Alpar et al. 1984 : Mendell 1991 : Andersson et al. 20061 ) . We should also consider the interaction between the 



vortices and the nucle i in the neutron star crust, the potential pinning of the vortices to the lattice (jPonati fc Pizzochero 



2006 



Avogadro et al. 20081) and th e extent to which the vortices exhib it creep in the presence of a rotational lag ( Anderson fc Itoh 
19751 : 1 Alpar et al.ll 19841 . Il989l : iLink et al.|[l99i : Ulpar et al.lll993l ) . 

In the following, we will focus on the role of the mutual friction and the glitch event itself. The vortex pinning will be dealt 
with in a very simplistic fashion. We will simply assume that the vortices are either perfectly pinned or completely unpinned. 
In such a model, a glitch would proceed as follows. Assume that the superfluid vortices form a uniform, straight, array aligned 
with the rotation axis (that is, we are not considering potential vortex tangles and turbulence I Andersson et al.l 2007 ) , whi ch 



would make the problem quite a lot harder since the model would have to contain local information (jPeralta et al 



20051 )) 



In this case vortex pinning simply fixes the number of vortices per unit area. This in turn fixes the neutron fluids angular 
momentum, so the superfluid component rotates at a constant rate. If we assume that the charged fluid is locked to the crust 
via magnetic effects then the vortices will be rotating with the charged fluid component. As the crust spins down due to 
the electromagnetic torque, a velocity difference will build up between the two constituents. This will lead to an increasing 
Magnus force acting on the vortices. Eventually, when some critical lag, AQc, is reached this force will be strong enough to 
overcome the nuclear pinning and the vortices are suddenly free to move. At this point the vortex mutual friction becomes 
relevant and serves to transfer angular momentum between the two components. This becomes the mechanism by which the 
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two components couple and the lag decays. The crust spins up leading to the observed glitch jump. If the system relaxes 
completely, the end state should be such that the two components rotate at the same rate. The glitch event itself is relatively 
sudden. The best resolved event to date is the so-called Christmas glitch in the Vela pulsar, where the glitch rise time was 
shorter than a few tens of seconds ( Dodson et al. 20021 ) . In other words, the angular momentum is transfered to the crust 



in less than a few hundred rotation periods. On a longer timescale one would expect the vortices to repin. After all, in the 
relaxed state the Magnus force is absent (or at least very small). The repinning should determine the long-term relaxation 
of the glitch, i.e. the spin evolution on timescales longer than tens of seconds. In order to model this phase one would likely 
need to account for vortex creep. Eventually, the system will reach a state where the rotational lag increases, and the pulsar 
may glitch again. 

We focus on the glitch event itself, i.e. the short term evolution following global vortex unpinning. During this phase one 
would expect the main dynamics to be determined by the mutual friction force. Hence, we do not have to con sider either the 



electromagnetic torque or the vortex creep. These are, of course, important in a complete glitch model, c.f. iLarson fc Link 



1 2OO2I) 'or an interesting discussion of the relaxation phase, but since they require dynamics on rather different timescales it 



is natural to (at least initially) consider the different phases separately. Our main aim is to compare a simple "global" glitch 
model to a numerical evolution that takes into account the actual hydrodynamics. 



4.1 The evolution equations 

To model a glitch event we need to account for the mutual friction force. As discussed by Andersson et al. ( 2006l l this means 
that we consider the evolution equations 

£t = f 4 + («r + exwD + V, (0 + Mx) + e^wYVrvi = ~ B^B' e^.u^iw';^ , (51) 

\ Ot J Px Px 

where B and B' are the mutual friction parameters ( Andersson et al. 20061 ). We have assumed that the vortex array remains 
straight and that the rotational lag remains orthogonal to the rotation axis. This is the simplest assumption. If we were to 
relax it, we would have to consider possible precession of the system and the mutual friction force would be more complicated. 
In the case that we are considering, we have 

ujI = 20^, + £n {20.^ - 2n'l) . (52) 
Generalising the prescription from the previous section, we can find the energy equation for each constituent. This leads 

to 

a^x 

dt 



PxU^ 



K? Pn I I yx 
. px ^ 



K,l Pn k I 

B — ejfc;aj,iTOyx 

Px 



dV. 



(53) 



The second term in the integral will vanish as vf is parallel to wf^. Written in terms of the moments of inertia I^i, cf. H26p. 
the total change in energy is given by 



^ = ^ <j ^In^n [fin + £„ (f^p - ^n)] + ^/pS^p [Op + £p (S^n - ^p)] } = -B \uJn\ (^p - Qn)' In < 



(54) 



That is, the mutual friction leads to a loss of kinetic energy (as expected). The equilibrium (minimum energy) state is reached 
when the two fluids are rotating together (Qp — Qn)- 

We can also calculate the global change in angular momentum. Focusing on the ^-component of the angular momentum 
we find 



Noting that 



and 



= j p^fiijkx^ £^ dV ^ j eijfcX^ l^iSpn |a;n| u)yx — S'pne''''"'*Jr«'m j dV . 



j klm n vx j n yx i n vx rx c 

€ijkX € OJi Wm = ^xf^z "^j ~ ^ OJj = (J lOr Z = Z , 



7 k 7 klm r^vx 

we can write the change in angular momentum in terms of the constituent moments of inertia to get 



^ijkX Wy^ — €ijjiX € iZ^ — "'j X Xj j ■ 



The total angular momentum is given by 



dXi__dJf_ dJl 

dt ^ dt ^ dt 



0. 



(55) 

(56) 
(57) 

(58) 
(59) 



Hence, the angular momentum is conserved. This is as expected since no external torques have been accounted for. It is worth 
noting that the B' coeflicient does not feature in the final equations. 
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4.2 An explicit solution 



We will now solve the global evolution equations for the system. To do this, it is useful to rewrite them in terms of the 
rotational lag and a quantity directly related to the total angular momentum. The conservation of angular momentum makes 
the latter variable trivial to deal with, while the lag is a key (more or less directly observable) quantity. 

From l|30p . (I59p and assuming that the moment of inertia of each constituent remains constant, we define V from the 
angular momentum such that 

/nf2n + /pf2p=/V = 0. (60) 

Then it is obvious that V remains constant during the evolution. From (|30|) and (|58[) we next find that the rate of change of 
the lag W = f2n — fip is given by 



(1 - £)W = Qn - fip 



where we have defined e ■ 



£n + £p- From (|60|l we can rewrite V as 
V 

Substituting this into (|61|) gives 



y n„ + = + ^ (a. - w) = - . 



(1 -e)W = -2B 



v + 



-w. 



(61) 



(62) 



(63) 



This is the stage at which the change of variables helps us. Because V is constant, equation (|63p is separable. Straightforward 
integration, assuming that the glitch occurs at time t — and defining Wo = W(0), leads to 



W = V 



V 

Wo' 



(="'-') 



The spin-up time r is given by 



(1 - e)/p 



(64) 



(65) 



2BIV 

For practical purposes it is better to express V in terms of the initial conditions. Defining the initial rotation of the protons 
as fio we easily find V at time t = 0. From (|60|) it follows that 

/p 



V : 



J (Jnfin + Ip^p) = y 



V 



^Wo 



This rearranges to give 



In 



Wo + Qo ~ Qo , 



which should hold since Wo ^ 1- We then arrive at the final result 



W ^ ^0 



^0 t/T 

Wo"" 



En 



Woe" 



-t/T 



where 



(1 - e)Ip 



2310.0 

It is easy to show that the observed component (the protons) evolves according to 



fio + y Wo ( 1 



(66) 
(67) 
(68) 
(69) 
(70) 



4.3 Matching Observational Data 

The model we have described is obviously quite simplistic. Most importantly, the assumption of constant parameters (which 
allowed us to carry out the integration over the body of the star in the first place) is quite unrealistic. Having said that, it 
would not be surprising if the final model were to retain some of the bulk dynamics of the more complex system. Of course, 
the various quantities in, for example, (|69p must be take to represent "body averages" in some sense. Moreover, the model is 
only relevant on the relatively short timescale of the glitch jump itself. In order to describe the subsequent long-term evolution 
we would need to include both the the magnetic spin-down torque and the repinning of the vortex lines to the crust. The 
latter could possibly be accounted for by "switching off" the mutual friction "gradually". That is, one could simply take B to 
be time-d e pende nt, refiecting the amount of vortex pinning or the nature of the vortex creep. This idea has been considered 



bv lSidervl (|2008l l. and we think it would be interesting to develop it further. 

Despite these caveats, it is interesting to consider how observations may constrain the various parameters. Let us consider 
the scenario where the observed component represents a small fraction of the total moment of inertia. This would be the case 
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for a typical neutron star crust coupled to a large superfluid reservoir in the core, when we may have Jp/J ~ 10~^ Then the 
observed glitch jump would be 



f2p — Oo 

no 



/n Wo 

/ fio 



Wo 
^0 



(71) 



That is, Wo would correspond (more or less directly) to the observed glitch size. At the same time, the available constraint 
on the glitch rise time can be compared to the spin-up time of the model. Let us, for simplicity, impose the constraint that 
the glitch happens in less than 100 rotations. Then we need 

(1 - e)h 



rflo 



2BI 



< 10^ 



We can rewrite the entrainment factor in terms of the effective proton mass in the usual way jPrix et aLlbood ). Then 



and we need 



or, for the suggested moment of inertia ratio. 



UlRj^ < 10^ 
mp 2BI 



e > 5 X 10 



-5 '"-p 



(72) 



(73) 



(74) 



(75) 



This constraint is not very severe. In particular, the canonical value B ~ 10 * ( Alpar et al.lll984 : Mendel! 1991 : Andersson et al 



2006ll lies within the required range. Howev er, if we use the contraint of a spin-up of the order of a day suggested by a par- 
tially resolved glitch in the CRAB pulsar jLvne et al.lll992l '). then the mutual friction parameter would be constrained to 
being weaker than B ~ 10"^". This result would not accord well with our current mutual friction models, likely illustrating 
our level of ignorance about the relevant physics. 

In principle, we would now want to consider t he case when the superfluid component represents only the free neutrons 
in the crust. That is, when we have In/ 1 ~ 10-2 jLink et ahlliiigl ). However, the model that we have developed does not 
immediately apply to this situation. This is obvious since we have assumed that two-fluid hydrodynamics applies throughout 
the system. In order to adress the crust superfluid problem we would have to add a component representing the single fluid 
core, and ensure that it is coupled to the two-fluid region in a suitable way. In particular, this core component would affect (|59p . 
This generalisation is complicated by the fact that we would have to add appropriate boundary conditions at the interfaces. 
As our main aim is to compare the averaged model to the detailed hydrodynamics we will leave consideration of models with 
distinct superfluid regions for future work. 



4.4 Energetics 

In the next section we will consider the hydrodynamics associated with a (core) glitch event. A key motivation for this 
discussion is the need to understand the actual details of how a macroscopic glitch is triggered (presumably through a large- 
scale instability) and how the system evolves once vortices become unpinned. By modelling the required hydrodynamics we 
hope to understand the nature of glitches better. We should, for example, be able to establish to what extent the simple "bulk 
model" we have discussed represents the behaviour of a true two-fluid system. We can also address other interesting questions, 
concerning for example the modes of oscillation that are excited in a glitch. This is an interesting question to ask because, 
first of all, there may be additional variability in the glitch event and, secondly, the fluid oscillations may be associated with 
gravitational-wave emission. It is obviously relevant to try to understand the nature of this gravitational-wave signal and 
estimate its amplitude. One should probably not expect glitching pulsar to be supreme gravitational-wave sources, but these 
estimates are nevertheless interesting. Most importantly, since glitches are common in young pulsars (and magnetars) it may 
be "reasonable" to assume that the corresponding level of energy is associated with regular dramatic events in a neutron 
star's life. 

Let us therefore consider the energetics of the problem. In past studies, it has been common to estimate the energy 
available for radiation based on a single component model. In that case the total kinetic energy and angular momentum are 
(obviously) given by 

E = i/n^ , and J ^m. (76) 

Assume that a glitch of size Afi results from a change in the moment of inertia A/. This would represent a "starquake" in 
an elastic star. Then, assuming that the total angular momentum is conserved, it is easy to show that the available energy is 

AEi^]^mAQ.. (77) 



As discussed by, for example, lAndersson fc Comeij l|200lf ) this estimate suggests that pulsar glitches may be of interest for 
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future generations of gravitational-wave astronomers. Ifowever, if we consider the two-component model we get a rather 
different picture. In this case, for constant 7x, the conservation of angular momentum in the glitch leads to 

Afin^-^Afip. (78) 

That is, the superfluid (neutrons) spin down as the crust (protons) spin up. Estimating the available energy, we find that 

ASa ~ i/p(A!:!)' . (79) 

Here, Af2 = Ailp and we have assumed that Jp ^ In- For typical parameters, I-p/I ~ O.f and Af2/f2 ~ fO~^, we see that 

A_B2 ~ 5 X 10"**A£i . (80) 

In other words, in the two-component model the energy available for radiation is much smaller than in the starquake case. 
Even though it is not clear how the estimate will change if we account for the energy radiated as heat, changes in internal 
and potential energy etcetera, it is clear that the result is rather pessimistic. If the estimate is taken seriously, and glitches 
really represent a transfer of angular momentum as in the two-fluid model, then the gravitational-wave signal from a pulsar 
glitch is unlikely to be detected by any future generation of detectors. Of course, our level of understanding of this problem 
is still rudimentary. We need to improve our models considerably if we are to make more reliable estimates. The simulations 
that we will now discuss provide an important step in this direction. 



5 HYDRODYNAMICS MODEL 



The "bulk" model that we have discussed so far is able to describe some key properties of glitches. However, the model has 
obvious restrictions. In particular, it does not provide any information whatsoever about the actual hydrodynamics of the 
event. By focusing on solid body motion we are obviously considering only the averaged dynamics. In principle, one would 
expect the result to be relevant when the dynamics is much slower than, say, the speed of sound in the fluid. However, it is 
clear that we need to move beyond the averaged model if we want to understand issues like the trigger mechanism for glitches, 
consider potential gravitational-wave signals etcetera. It is thus natural to consider the hy drodynamical aspects of the glitch 
proble m. T o do this, we have extende d the numerical code that w as recently developed bv lPassamonti. Haskell Andersson 
l|2009h [see IPeraha et al] (|2005l . bood ): IPeralta fc MelatosI l|2009h for a parallel effort]. Within the two-fluid framework, we 
evolve perturbations of rotating, superfluid Newtonian stars in time. The new version of the code includes the effects of mutual 
friction and the perturbed gravitational potential (i.e. we are not working in the Cowling approximation). We initiate the time 
evolution with suitable conditions that mimic a pre-glitch configuration, where the two fluids rotate uniformly with different 
velocities. Assuming that the vortex pinning that is required to reach this state is instantaneously broken, we can evolve the 
system. This allows us to determine the mutual friction damping, extract the associated gravitational signal and infer the 
oscillation modes that are excited during a "glitch". In particular, we can test the analytical formula for the global spin-up 
timescale r, equation (|69p . 

The main motivation for our perturbative treatment is to consider stellar models where the relative velocity lag between 
neutrons and protons is very small as a deviation from stationary equilibrium configurations. These background models, which 
are such that the two fiuids co-rotate, are in /3-equilibri um and coexist throughout the star's volume, can be constructed by 
extending the standard self-con sistent field method of Hachisu (|l986 ). The deta i ls of this method, and its application to 
superfiuid stars can be found in Yoshida fc Eriguchi (2004) and Passamonti et al. ( 20091 ) . In our current models the crust is 
neglected (for simplicity). Our aim is to continue to add key physics to the model, and we plan to consider crustal effects in 
future work. 



Non-corotating configurations can be determined using the perturbative approach developed by lYoshida fc Eriguchi 



(2004), where the deviations from co-rotation are numerically computed by means of a variation of the self-consistent field 
method. We extend this numerical approach to study different superfluid equations of state and generate rotating stellar 
sequences with constant mass. These non-corotating deviations are then implemented in the hydrodynamical code as initial 
conditions and evolved in time with a system of perturbation equations. This system is formed by the linearised versions of 
the two-fluid mass conservation equations, the two Euler-type equations ()5ip and the Poisson equation for the gravitational 
potential. Technical details on the constr uction of the initial data and the time domain numerical code will be provided 
elsewhere ( Passamonti fc AnderssonI 20091 ) . Here, we report only results that can be directly compared with the analytic 
expressions from the previous sections. 

In a two-fluid model, the relative motion between protons and neutrons can be approximately decomposed in co- and 
counter-moving components. With an appropriate choice of perturbation variables, w e can study the effects of these two degrees 
of fre edom on the oscillation spectrum and the factors that generate their coupling ( Andersson et al. 20091 : Passamonti et al. 
20091 ). 

The perturbation equations must be completed with an equation of state (EoS), which can be described by the energy 
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Table 1. This table provides tlie main parameters of the corotating background models for both the A and C sequences. The first column 
labels each model, while the second and third columns give, respectively, the ratio of polar to equatorial axes and the angular velocity 
of the star. In the fourth column, the rotation rate is compared to the Kepler velocity Qj^ that represents the mass shedding limit. The 
ratio between the rotational kinetic energy and gravitational potential energy T/|iy| is given in the fifth column. In the sixth and seventh 
columns we show the moment of inertia of the proton and neutron fluids, respectively, while in the eighth column we provide the stellar 
mass. All quantities are given in dimensionless units, where G is the gravitational constant, po represents the central mass density and 
Req is the equatorial radius. 



Model 


Rp/ Req 






T/\W\ X 10^ 


Ip/iPoRlq) 


Ir./{poRlq) 


M/{poBlq) 


AO 


1.00000 


0.00000 


0.00000 


0.00000 


0.03328 


0.29951 


1.2732 


Al 


0.99792 


0.05913 


0.08156 


0.05802 


0.03319 


0.29868 


1.2701 


A2 


0.98333 


0.16675 


0.22999 


0.38482 


0.03253 


0.29278 


1.2479 


A3 


0.95000 


0.28799 


0.39627 


1.16918 


0.03102 


0.27915 


1.1967 


A4 


0.93333 


0.33081 


0.45629 


1.56885 


0.03025 


0.27224 


1.1709 


A5 


0.90000 


0.40268 


0.55543 


2.38295 


0.02869 


0.25822 


1.1186 


CO 


1.00000 


0.00000 


0.00000 


0.00000 


0.01906 


0.24657 


1.0826 


CI 


0.99792 


0.05586 


0.08403 


0.04561 


0.01900 


0.24584 


1.0798 


C2 


0.98333 


0.15764 


0.23716 


0.36682 


0.01858 


0.24064 


1.0601 


C3 


0.95000 


0.27145 


0.40837 


1.11303 


0.01760 


0.22861 


1.0146 


C4 


0.93333 


0.31246 


0.47006 


1.49236 


0.01711 


0.22251 


0.9915 


C5 


0.90000 


0.38006 


0.57177 


2.26285 


0.01610 


0.21009 


0.9447 



functional (|16p : 

E = E {pn,pp,wlp) , (81) 
where we have replaced the number densities rix with the mass densities px (assuming for simplicity that the neutron and 



proton masses are equal, Trip = mn). When the relative velo city between the two fluids is small, equation (|8ip can be expanded 
in a Taylor series (|Prix et al-lbooi IPassamonti et ahlbood ) . Then we have 

E = Eo {pn, Pp) + ao (Pn, pp) wlp + O (Wnp) , (82) 
where the entrainment function, ao, and the bulk equation of state, Eo, can be independently specified on the corotating 
ba ckground, where = 0. In this paper, we consider two sets of polytropic superfluid equations of state already used 
by IPassamonti et al.l 1 20091 ) . Despite their simplicity, these models are very useful for investigating the effects of the different 
physical parameters on the oscillation dynamics. In particular, one of these two classes of equation of state describes models 
with composition gradients that couple the co- and counter-moving degrees of freedom. 
The first equation of state is given by the following expression: 



Eo 



K 



-Pl 



2Ka 



-pnPp + 



K[l + a 



(l + 2cr)a::p] 2 
-Pp , 



(83) 



1 - (1 +t7)a;p''" 1- (l + o-)a;p'""^^ ' 3: p [1 - (1 + a ) j 

where K and the proton fraction Xp are taken to be constant. As discussed bv lPrix et al. (2002), the parameter a is related to 
the "symmetry energy" of the EoS. For this model, we have constructed a sequence of co-rotating axisymmetric config urations 



which do not have composition gradients. These non-stratified stars correspond to the A models used by Passamonti et al 



( 20091 ). Table [U gives the main quantities of some rotati ng models that we consider in this paper. Details for more rapidly 

rotating models can be found in Passamonti et al. ( 20091 ) . 

The second analytical equation of state, t hat describes stratified superfluid stars, can be written I Prix fc Rieutord 20021 : 
Andersson et al. 20021 : Passamonti et al. 20091 ): 

Eo = kn pT + fcp p7 ■ (84) 
In this paper we determine a sequence of rotating stars, whose non-rotating member is Model III of 'P rix fc RieutordI (|2002l ). 
For our polytropic models we use 7n = 1.9 and 7p = 1.7, while the coefficients fcx are given, in units GReqpl"'^ by K = 0.682 
and ftp = 3.419, respectively. Recall that G, R^q and po are the gravitational constant, the equatorial radius and the central 
mas s density , respectively. Imposing /3-oquilibrium on the corotating background model, we can determine the proton fraction 
as (|Prix fc Rieutordi200a : iPassamonti et al., .2009'): 



1 + 



(7pfcp) 

(7iifcn) 



-JV„ 



(85) 



For this modelj_a se- 



where A^n and A^p are the neutron and proton polytropic indices of the EoS, defin ed as iVx = (7x — 1) ^ ■ . 
quen ce of rotating stars can then be determined via the self-consistent field method ( Passamonti et al.|2009l: Passamonti fc Andersson 
20091 ). We refer to these stars as models C, in order to distinguish them from the models B used by IPassamonti et al.l (|200^ 



See Table [T] for some of the rotating configurations for this model. 
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5.1 Glitch initial data 

According to the two-fluid model, a glitch brings a star with an initial velocity lag between protons and neutrons to a co- 
rotating configuration. The observed glitch jump is very small, Ail/fi <C 1, and can be associated with angular momentum 
transfer from the superfluid neutrons to the proton component. Regardless of the physical mechanism that generates the 
original velocity difference, we can describe the initial conditions for a glitch as an axisymmetric configuration where protons 
and neutrons rotate with a small relative velocity. If we assume tha t this relative rotation is al igned with the "background" 



rotation axis, we can use the perturbative approach developed by lYoshida fc Eriguchil (|2004l ) to determine the difference 
between the initial non-corotating configuration and the final corotating background. 

As we have already discussed, the observed variation of the star's angular velocity during a glitch can be associated 
with the proton velocity. In fact, due to the magnetic field coupling between the crust and core protons, we can assume that 
all charged particles corotate. Meanwhile, the velocity of the superfluid neutrons must be determined from the conservation 
laws. If we assume that the angular momentum of the two-fluid system is conserved during a glitch, then the initial relative 
amplitude between the proton and neutron perturbations can be determined from 

^ A Jx = , (86) 

X 

where the perturbed angular momentum of the x component is given by 

AJx = /xA!^x + AJx!^, (87) 

and 7x is the moment of inertia of the x fluid in the corotating conflguration, which rotates with angular velocity Q.. We 
determine the perturbation A7x from 

A7x= / (5px(r'sin6»)^dr', (88) 

where (5px is the Eulerian perturbation of the mass density. From equation (|86p we then have 

JnAJ^n + JpAfip + (A7„ + Alp) n = , (89) 
which leads to the following expression for the initial fluid angular velocities: 



(90) 



Afip _ AO. 

AQn = (JpA^p + A/f7) , (91) 

^11 

where A7 = Alp + AI^. 

For any corotating background we can determine t he non-corotating correction s to the mass density 5p^, the chemical 



poten tial Sjl^ and the gravitational potential 5<j) via the lYoshida fc Eriguchil (|200J) approach, see iPassamonti fc Andersson 



:br more details. In a spherical coordinate basis, the initial velocity fields are then given by 

<5wx = -Af7x ■ (92) 

This initial data is close to that considered in the bulk dynamics model of the previous section. Hence, we can meaningfully 
compare the results of our time evolutions to the averaged results. This comparison will give us a better idea of the validity 
of the solid-body approach. Of course, by solving the hydrodynamics problem we will be able to proceed beyond the averaged 
model and discuss other interesting issues. 



5.2 Spin-up time 

Let us flrst compare the mutual friction damping extracted from the time-evolution to the body averaged analytical for- 
mula (|69|l . Considering the sequence of non-stratifled models A, we test the dependence of r on the mutual friction strength 
B, the background angular velocity Q, the moment of inertia ratio Ip/I and the entrainment parameter e. For the corotating 
background A models given in Table [U we determine the initial condition for the time evolution code that corresponds to a 
"typical" glitch jump AQ/Q, = 10~^. The related neutron velocity lag is then given by equation (|9ip . It is worth noting that 
the method used to construct the initial data is linear in the two fluid velocities. This means that the results can be rescaled 
to other values of the glitch size. 

If we exclude the mutual friction term, the numerical evolutions preserve the initial lag between protons and neutrons. 
At the same time, the initial conditions excite low level oscillations. This is expected as we are mapping a non-corotating 
axisymmetric conflguration onto an axisymmetric corotating background. When mutual friction is switched on, the counter- 
moving motion is damped, leading to a corotating conflguration on a timescale of order r. We determine the spin-up timescale 
by monitoring the ip component of the velocity difference Wp^, assuming that it depends on time as iOp„ = Wp^ (t = 0) e~''^^. 
By taking the natural logarithm, we can extract r from a linear flt of the time evolved data. 
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Figure 1. This figure shows the spin-up time t as a function of the inverse of the mutual friction parameter B (left panel) and the 
background rotation rate of the star il (right panel) for the sequence of rotating models A. The solid lines show the behaviour predicted 
by equation 16911 . while the symbols (see legend) represent the values of the spin-up time extracted from the hydrodynamical simulations. 
The physical quantities are given in dimensionless units by using the gravitational constant G and the central mass density po- The axes 
use logarithmic scales. All the five models A1-A5 shown in this figure have both vanishing symmetry energy term a and entrainment 
parameter e. In the left panel, the proton fraction is fixed to Xp = 0.1 and the mutual friction is varied. Meanwhile, in the right panel we 
show three sequences of rotating stars with the same mutual friction strength S = 0.1, but with three different values of proton fraction, 
namely Xp = 0.1, 0.05 and Xp = 0.01. In all cases, the numerical values of spin-up time show a good agreement with the analytical result. 



In Figs. [l]and[2l we show the agreement of the numerical results with the analytical formula (|69|l . In the left panel of 
Fig. [1] we consider several rotating A models with vanishing symmetry energy and entrainment and with constant proton 
fraction Xp — Ip/I — 0.1. The models have different mutual friction strength, controlled by the parameter B that we take 
to be constant throughout the star. The expected linear dependence of r on the inverse of the mutual friction parameter is 
confirmed by the hydrodynamical simulations. In the right panel of Fig.[T] we fix instead the mutual friction strength, B = 0.1, 
a rather large value, and study three sequences of rotating models with different proton fraction, respectively Xp = 0.01, 0.05 
and 0.1. The mutual friction damping time exhibits the expected linear dependence on fi"^. For models Al and A2 we test 
also the dependence of r on the entrainment parameter e, see Fig[2] In this case, the other stellar parameters are, respectively, 
Xp — 0.1, B = 0.1 and a — 0. The linear dependence on 1 — e is clearly confirmed by the evolutions. According to equation (|69p . 
the damping time should not depend (explicitly) on the symmetry energy, a. We have carried out simulations with different 
a to confirm this result. 

Next, we consider the stratified models C. The aim is to establish to what extent equation (|69p still provides accurate 
results for the spin-up time. On the one hand, one may not expect this to be the case since the various parameters in the 
model are no longer uniform. On the other hand, the simple prescription could still work provided that the parameters are 
interpreted in a body-averaged sense. We consider initial configurations with Aft/Q = 10~^ and determine the non-corotating 
corrections using the method discussed in Section FS.ll For the three rotating models C1-C3, the dependence of the numerical 
damping time on the mutual friction parameter B is shown in the right panel of Fig. (2] Again, the results agree well with 
the analytical formula (|69|) . However, a closer examination of the data reveals that equation (|69|) is more accurate for the 
non-stratified A models. For the first three rotating models of the two sequences A and C, we show in Fig. [3] the relative 
deviation between the numerical and analytical spin-up time for different mutual friction strengths. Comparing models with 
the same axis ratio, the error is generally smaller for models A (filled symbols). 

In general we find that the agreement between the numerical and analytical spin-up time is better for slowly rotating mod- 
els that have strong mutual friction. In these cases, the damping of Wp^ is less contaminated by the excitation of axisymmetric 
oscillations, and the numerical extraction of r is more accurate. 



5.3 Gravitational Waves 

So far we have discussed "global" dynamics, e.g. how the two components in the system relax to a co-rotating configuration 
due to the mutual friction. We now turn to the actual hydrodynamics, and consider to what extent this kind of glitch event 
radiates gravitational waves. We have already established that these events are unlikely to be strong emitters of gravitational 
radiation. However, it seems inevitable that they should radiate at some level and it is important to establish what that level 
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l-l l/B 

Figure 2. In this figure we compare tlie numerical spin-up time r with the analytical formula I I69II . The values extracted from the 
numerical code are shown as open symbols (see legend), while the solid lines denote the analytical t. In the left panel, we consider the 
two models Al and A2, with proton fraction Xp = 0.1, vanishing symmetry energy ct = and constant mutual friction parameter B = 0.1. 
The dependence of the numerically determined r on the entrainment parameter agrees very well with the analytical result. In the right 
panel, we show (on a log-log scale) the damping time t as a function of for the three models C1-C3 with vanishing entrainment. 
The agreement between the numerical and analytical spin-up times is still good, although less accurate than for the A models. See Fig. |3] 
and the main text for further discussion. 



may be. In particular, since there are a number of glitching pulsars in the Galaxy. The energy involved in glitches indicates 
how energetic "typical" events in a mature neutron star may be. In that way, these events provide interesting benchmarks 
for gravitational-wave modelling. It is, however, important to state already from the outset that we are not expecting the 
mechanism that we are considering here to lead to detectable signals. This is essentially because of the high level of symmetry 
in the initial and final configurations, and the fact that we are assuming global vortex unpinning. The result may be quite 
different if we were to model localized unpinning events. Although of great interest, this problem is unfortunately beyond the 
reach of our current computational technology. 

We will focus on the gravitational signal associated with the I — 2 axisymmetric oscillations that are excited in our glitch 
evolutions. At the linear perturbation level, the initial data excites a number of the neutron star's oscillation modes. Hence, 
a key question concerns which modes we expect to be present in the gravitational signal. For a single fluid star, the general 
mode classification is based on the main restoring force that acts on the displaced fluid elements (Cowling 1941). In this work 
we consider "acoustic modes" that are mainly restored by pressure variations. Since any perturbation of a spherical star can 
be decomposed in vector harmonics, an oscillation mode can be labeled by the harmonic indices (/,m) associated with the 
spherical harmonics Yi^{6, <j)). This description can be extended to rotating stellar models, as the modes can be tracked back 
to the non- rotating limit. For any value of {l,m), the oscillation modes can be ordered by the number of radial nodes in their 
eigenfunctions. For acoustic modes, we then have the fundamental mode 'f with no nodes and the series of pressure modes 
'pj with i nodes. 

Our glitch model leads to axisymmetric initial data. Therefore, we can only excite the family of axisymmetric modes, 
with m = Q. The quadrupole, I — 2, oscillations are expected to be dominant in the gravitational signal. However, in rotating 
stars the gravitational-wave spectrum can also contain I — "quasi-radial" oscillation modes. In the non-rotating limit, 
these become purely radial modes, which do not generate gravitational radiation. The quasi-radial fundamental mode will be 
denoted by F and its i overtones by H^. 

In a two-fluid model the oscillation spectrum is richer, as the displaced fluid elements can now oscillate in phase and 
counter-phase. The comoving degrees of freedom generate oscillation modes that are similar to those of a single fluid star, 
and are referred to as "ordinary modes". The counter- moving degree of freedom produce a new class of modes known as 
"superfluid modes". In our discussion, ordinary and superfluid modes will be labeled with an upper index, for instance the 
/ = 2 fundamental ordinary mode will be expressed as ^f°, while the superfluid mode as ^f 



In a two-fluid star, the gravitational radiation is entirely generated by the co-moving degree of freedom (jAndersson et al 

20091 ). We determine the gravitational strain from the standard quadrupole formula (|Thorneiil980i ): 



hf = ^i±l T,f'^\ (93) 
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where the {I, m) 



(2, 0) pure spin tensor harmonic is given by 

1 /15 



8 



The mass quadrupole moment for the two-fluid star is defined by l Andersson et al. 20091 ) 



15 

where P^" is the associated Legendre polynomial: 



Scos-'e- 1 



(94) 



(95) 



(96) 



In equation (|93|1 . the numerical calculation of the second order time derivative of the quadrupole moment may lead to in- 
accuracies. However, the gravitational-wav e extraction can be improved by using the dynamical equations and replacing 
the time derivative by spatial derivatives ( Finn fc Evans 19901 ). In our perturbative analysis, we have found accurate re- 
sul ts already with the mom e ntum -formula, where only first time derivatives appear. More details and tests will be given 
by Passamonti fc AnderssonI ( 2009l l . 

We rewrite equation (|93[l as follows: 
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where the quantity At^ is given by (jPassamonti fc Anderssonll2009 
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The energy radiated in gravitational waves can be determined from the following equation (|Thornelll980l ') 
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and A^'^ is its Fourier transform. The energy spectrum is then given by 
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In our analysis, we focus on the two rapidly rotating models A2 and C2, which rotate at a significant fraction of the mass 
shedding limit ^1/Q,k ~ 0.23 (see Table [T] for more details). For a typical neutron star with mass Af = IAMq and radius 



Re 



10 km, the rotation period of the A2 and C2 models is about P ~ 3 ms. These models are therefore rotating much 



faster than all known glitching pulsars. We will demonstrate that the generation of gravitational waves is very small even 
for these relatively rapidly rotating models. The result can be considered as an upper limit for real glitching neutron stars. 
Moreover, we will show that the gravitational signal of slower rotating models can be determined by a simple rescaling of the 
A2 and C2 signals. 

At the linear perturbation level, the entrainment parameter e can be chosen independentl y from the bac kground model. 
Recent work suggests that this parameter can assume values in the range 0.2 < e < 0.8 (IChame]l2008l). The effect of 



the en trainment on the oscillation spectrum has been extensively studied by Prix fc RieutordI (|2Q02l ) and 

( 20091 ). The oscillation modes that are mainly affected by e are the so-called "superfluid" modes, which are associated with the 
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counter-moving degree of freedom. Simple relations between the frequ e ncies o f the superfiuid acoustic a nd inertial modes an d 
the entrainment parameter e have been obtained by Andersson et al. ( 20091 ): Passamonti et al. ( 20091 ): Haskell et al. ( 20091 ) . 
In this work we show results only for the e = 0.5 case, as apart from the spectral properties discussed above, we did not find 
any qualitative difference in simulations using other values of this parameter. 

With the method discussed in Section 15.11 we set up initial data that mimic the configuration of an axisymmetric 
glitch. The initial value for the proton and neutron angular velocity can be determined from equations (I90|l - (|91|l . once we 
fix the glitch size Ail/Q and solve the stationary equations for the background model. We will consider the case of a large 
glitch, where AQp/Q = 10~^. This means that, due to angular momentum conservation, the neutron fluid slows down with 
AQn/^ = -1.11 X 10"^ for model A2 and AQn/^ = -7.74 x 10"** for model C2. Note that we use a flrst order perturbative 
framework, where for a given corotating background model the results of the time evolutions are linear with respect to the 
parameter AQp/fl. Hence, the gravitational-wave strain can be rescaled to any desired glitch magnitude. Furthermore, for 
slower rotating models that have the same glitch size, Afip/f2, we expect the perturbations and the gravitational strain to 
exhibit a quadratic dependence on the background rotation rate fl. Our numerical simulations reproduce this behaviour when 
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Figure 3. In this figure, we estimate tlie agreement between the analytical spin-up time T from equation II69I I and the values extracted 
from the hydrodynamical code. We show the relative deviation At/t as a function of the inverse of the mutual friction parameter B. The 
horizontal axis is on logarithmic scale. The results for models A and C are shown as filled and open symbols, respectively (see legend for 
details). 



the stars have small rotational deformations, as in the case of the A1-A2 and C1-C2 models. Already for models A3 and 
C3, this scaling with is less clear and obviously it is not expected to hold for more rapidly rotating stars. In conclusion, 
from the evolutions of the A2 and C2 models, we can easily estimate the gravitational strain emitted by other slowly rotating 
models with different glitch size and background rotation. 

In Fig. [4] we show the characteristic strain he for the A2 and C2 models with Aflp/Q = 10~®. The results refer to a 
star with mass M = IAMq, radius Req = 10 km and with a low level of mutual friction. We consider an evolution that 
lasts for ~ 27.25 ms and extract the signal at 1 kpc0. This is roughly the distance to the Vela pulsar. The first key result 
in Fig. [4] is that the gravitational signal of the C2 model is about ten orders of magnitude larger than that of the A2 model. 
This is an enormous difference, given that these ought to be the same kind of events. The difference is due to the presence 
of composition gradients in the C models. In a stratified model, the co- and counter-moving degrees of freedom are coupled 
during the evolution. This coupling is crucial, since only the co-moving motion generates gravitational radiation. The initial 
data for the pre-glitch lag between neutrons and protons, in accordance with equations (|90p - (|91[) . represent a counter flow. 
Hence, in the non-stratified A models these conditions generate a purely counter-moving motion between the two components 
that does not produce any gravitational signal at all. The fact that the strain of model A2 is not completely zero in the left 
panel of Fig.|4]is likely due to numerical errors. In fact, we have established that the level of radiation decreases with increased 
resolution. The result shown for model A2 corresponds to an initial non-corotating configuration where the variation of the 
total rotational kinetic energy is AErot/Erot ^ If we increase the precision of the self-consistent field method that 

provides the initial data we can lower this value and consequently the gravitational signal converges to zero. Moreover, the 
numerical noise in the simulation excites some oscillation modes that are related to the co-moving motion. In the left panel of 
Fig.|4]we identify the fundamental I = 2 mode ^f, the first two pressure modes ^pi and ^p2, and the quasi-radial fundamental 
mode F with its first overtone Hi. 

Let us contrast the results for the non-stratified A2 model with the results for model C2. In this case, the chemical 
coupling between the two fiuids introduces a co-moving motion already in the initial data. The evolutions then generate a 
larger gravitational strain and several oscillations modes, like the fundamental I = and I = 2 modes and their respective 
overtones. In particular, in the right panel of Fig |4] we note that both ordinary and superfluid modes are present in the 
gravitational radiation. This is due to the coupling of the degrees of freedom (the oscillation modes are no lo nger purely co- or 



counter-moving). The mode frequencies of the non-ro tating model CO have been com pared with the results of lPrix fc Rieutord 



I 2OO2I ). The two results agree to better than 1.4% (see lPassamonti fc Anderssonlboool ). We have identified the oscillation modes 



It is worth noting that we have evolved the system for approximately ten rotation periods. We did not extend the evolutions because, 
in reality, one would expect the coupling of the two components to start playing a role at this stage. If this were not the case and the 
oscillation prevailed for the entire timescale of gravitational-wave damping, then the f-mode may last a few seconds. This would be about 
a factor of 100 longer than our evolution, meaning that the effective gravitational-wave strain could increase by perhaps an order of 
magnitude. It would still bo too weak to be detectable. 
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Figure 4. This figure displays the gravitational-wave signal generated by our hydrodynamical glitch simulations for the two models A2 
(left panel) and C2 (right panel). On the horizontal and vertical axes we plot the oscillation frequencies and the characteristic strain 
extracted at a distance of 1 kpc from the source. We consider a neutron star with typical mass M = 1.4Mq and radius Rgq = 10 km. 
Models A2 and C2 then correspond to stars with rotation period P = 3.07 ms and P = 3.00 ms, respectively. The other stellar 
parameters are e = 0.5 for the entrainment and a = for the symmetry energy. In the case of model A2 the proton fraction is constant, 
Xp = 0.1, while model C2 is stratified with central proton fraction Xp(0) = 0.1. The initial configuration corresponds to a large glitch 
with AClp/Cl = 10~^, as described in the main text. We run the simulation for about 27.25 ms, i.e. about 9 rotation periods, and 
neglect the mutual friction force. From the displayed results, the strong effects of the stratification on both the oscillation spectrum and 
gravitational- wave amplitude are evident. 



of the C2 model by carrying out simulations with difTerent values for the entrainment e and tracking t he superfluid mo des as 
the parameter changes. To this end, we have also used the analytical formulae determined bv iPassamonti et al. ( 20091 ) . 



6 CONCLUDING REMARKS 

We have discussed the dynamics of pulsar glitch events from two, complementary, points of view. First we constructed a simple 
model based on global "averaging" of the standard two-fluid equations including the mutual friction due to superfluid vortices. 
This analysis provides a more detailed derivation of the phenomenological relations that have been used in many discussions 
of glitches. In particular, our final relations clarify how the spin-up time depends on key parameters like the entrainment. 
The derivation also highlights the various assumptions and the restricted validity of the model. Anyway, for typical values 
of the parameters (see Sec. 14. 3|) . our model has a glitch rise time shorter than the upper bound set by current observations. 
The model provides a useful description of the actual glitch event, but it does not account for the subsequent long-term 
relaxation (on a timescale of days to months) of the system. A key conclusion from our discussion is that the late stages of 
evolution requires additional assumptions, most likely, concerning the repinning of vortices. Understanding this phase better, 
e.g. connecting it to the two-fiuid hydrodynamics and the averaged forces that act on the vortices, is an important challenge for 
the future. It seems clea r that vortex creep will play a central role ( Anderson fc Itoh 19751 : Alpar et al. 19841 . 19891 : Link et al. 



19931 : lAlpar et al.lll993l ). but this mechanism has not yet been discussed in terms of the macroscopic hydrodynamics. This 
issue needs to be addressed if we are to develop more detailed models of glitch dynamics. We definitely need to move beyond 
phenomenology. 

As a first step toward s hydrodynamic glitch modelling, we have extended the recent linear perturbation evolution code 
of Passamonti et al. I (|2009l ) to include the mutual friction and the perturbed gravitational potential. Initiated with pertur- 
bations that represent two fluids rotating uniformly at different rates, the numerical code shows how the system relaxes to 
co-rotation. We have analysed this relaxation in detail and demonstrated that the behaviour is accurately described by the 
phenomenological model, at least for non-stratifled stellar models. When the star is stratified (e.g. has varying composition) 
the relaxation deviates from the simple model. This is as expected, since the global model was derived under the assumption 
of uniform parameters. Of course, the numerical evolutions provide us with a useful tool for studying the behaviour of more 
complex stellar models. In addition, our time evolutions provide a first insight into the excitation of neutron star oscillations 
by glitches. Our results show that a set of axisymmetric modes are excited by the glitch initial data. These modes will radiate 
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gravitational wave^, and it is important to establish if the associated signals may be observable with future detectors. In this 
respect, our results are quite pessimistic. In the cases that we have considered, the gravitational-wave signal is too weak to be 
detectable (even with a third generation of detectors However, it is not clear that this is the final say on the matter. One 
should keep in mind that the gravitational-wave strain differs enormously for our two model configurations. The non-stratified 
model does not (in principle) radiate at all, while the stratified model leads to a qualitatively interesting (albeit weak) signal. 
The enormous difference between these results shows that we need to continue to refine our modelling. We obviously have 
to account for the variation of composition throughout the star, and consider the fact that superfluid components will only 
be present in specific density regions. We also need to understand the nature of the vortex pinning better. In our models we 
have assumed that the vortices unpin in a catastrophic global event. It is far from clear that this is the case in a real system. 
It could, for example, be that the unpinning is localized. This would make the event less symmetric which may enhance the 
gravitational-wave sign al. We clearly need to understand t he actual mechanism that triggers the glitches better. The superfluid 
instability discussed bvl ciampedakis fc Andersson ( 20091 ) is interesting in this respect, but we need to study this mechanism 
in more detail to establish to what extent it can operate in a real neutron star. 

To make progress we need to overcome a number of challenges. Yet, recent developments have provided us with interesting 
insights and (most importantly) computational technology that should allow us to study much more realistic neutron star 
models in the not too distant future. 
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